Goal: append IBAs to preliminary CA-BBA blocks, visualize, and calculate basic summary statistics.
library(dplyr)
library(here)
library(mapview)
library(tigris)
library(sf)
mapviewOptions(basemaps.color.shuffle = FALSE)
ca.blocks.clean <- readRDS(here("Data", "CABlocks", "CABlocks2024.RDS")) %>%
filter(WATER_SLIVER == 0 & CA_SLIVER == 0)
cntr_crds <- c(mean(st_coordinates(ca.blocks.clean)[, 1]),
mean(st_coordinates(ca.blocks.clean)[, 2]))Important Bird Areas from Audubon California were mapped onto atlas blocks. Data from here.
Multiple IBAs sometimes overlapped a single block. In these cases, we recorded the ID of all IBAs, but for visualization purposes we prioritized the most-localized IBA (the one with the smallest total area).
# Load IBA data
IBAs <- st_read(here("Data", "IBAs", "iba_polygons_public", "iba_polygons_public.shp")) %>%
st_transform(st_crs(ca.blocks.clean)) %>%
st_make_valid() %>%
mutate(IBA_area = st_area(.) %>% units::drop_units())
# Intersect with CA-BBA blocks, count number of IBAs by block, and record thier IDs
iba_intersect <- st_intersection(ca.blocks.clean, IBAs)
iba_intersect <- iba_intersect %>%
mutate(intersection_area = st_area(geometry)) %>%
as.data.frame() %>%
st_drop_geometry() %>%
select(BLOCKCODE, site_id, intersection_area) %>%
mutate(intersection_area = units::drop_units(intersection_area)) %>%
group_by(BLOCKCODE, site_id) %>%
summarize(intersection_area = sum(intersection_area), .groups = 'drop') %>%
group_by(BLOCKCODE) %>%
summarize(
IBA_COUNT = length(unique(site_id)),
IBA_ALL_ID = paste(unique(site_id), collapse = ", "),
IBA_MOST_COVERAGE = site_id[which.max(intersection_area)],
.groups = 'drop'
)
# For blocks with multiple IBAs, identify which IBA is the most localized (smallest total area). Also pull names of all IBAs in each block using IBA ids.
iba_intersect$IBA_MOST_LOCALIZED_ID <- NA
iba_intersect$IBA_ALL_NAMES <- NA
for(i in 1:nrow(iba_intersect)){
rel_ibas <- IBAs %>% st_drop_geometry() %>% filter(site_id %in% as.numeric(unlist(strsplit(iba_intersect$IBA_ALL_ID[i],split=", ",fixed=TRUE))))
iba_intersect[i, "IBA_MOST_LOCALIZED_ID"] <- rel_ibas %>% slice_min(IBA_area) %>% pull(site_id)
iba_intersect[i, "IBA_ALL_NAMES"] <- paste(rel_ibas %>% pull(site_name), collapse = "; ")
}
# Join IBA info back to full blocks
ca.blocks.clean <- ca.blocks.clean %>%
left_join(iba_intersect, by = "BLOCKCODE") %>%
left_join(IBAs %>% st_drop_geometry() %>% select(site_id, site_name, priority), by = c("IBA_MOST_LOCALIZED_ID" = "site_id")) %>%
rename("IBA_MOST_LOCALIZED_NAME" = "site_name",
"IBA_MOST_LOCALIZED_PRIORITY" = "priority")Blocks in IBAs: 4824 of 16392( 29.4%)